m 
o 

(N 



•4— > 



o 
o 



> 
o 

00 
00 

o 



X 



Dynamical phase transitions in supercooled liquids: interpreting 
measurements of dynamical activity 

Christopher J. Fullerton 1 and Robert L. Jack 1 

Department of Physics, University of Bath, Bath, BA2 7 AY 

We study dynamical phase transitions in a model supercooled liquid. These transitions occur in ensembles of 
trajectories that are biased towards low (or high) dynamical activity. We compare two different measures of 
activity that were introduced in recent papers and we hnd that they are anti-correlated with each other. To 
interpret this result, we show that the two measures couple to motion on different length and time scales. We 
find that 'inactive' states with very slow structural relaxation nevertheless have increased molecular motion 
on short scales. We discuss these results in terms of the potential energy landscape of the system and in 
terms of the liquid structure in active/inactive states. 



I. INTRODUCTION 

As liquids are cooled towards their glass transitions, 
their relaxation times increase dramatically, and the mo- 
tion of their constituent particles becomes increasingly 
co-operative and heterogeneous 1-3 . There are several 
competing theories that aim to describe these phenom- 
ena 4-7 , but neither simulation nor experimental data 
have so far proven sufficient to establish which (if any) 
can fully describe the supercooled liquid state. Recently, 
novel dynamical phase transitions have been discovered 
in glassy systems 8-11 : these are new results that can be 
used to test existing theories. These phase transitions 
take place in ensembles of trajectories (sometimes called 
s-ensembles), where the dynamical evolution of the glassy 
systems is biased towards low-activity states 9,12-14 . Since 
these phase transitions are dynamical in nature, they fit 
naturally with theories of the glass transition where dy- 
namical motion takes a central role 4,15,16 , but they can 
also be interpreted in terms of random first order transi- 
tion theory 5 , and are linked with properties of the energy 
landscape and its normal modes 7,17-19 . 

In this article, we discuss these dynamical phase tran- 
sitions and their associated ensembles of trajectories. 
We are motivated primarily by two previous studies 20 ' 21 
which provided evidence for such transitions in a model 
glass-former, composed of Lennard- Jones particles 22,23 . 
In the first study, Hedges et al. 20 measured the activity 
in this model through the mean square displacement of 
its particles. Biasing the dynamics with respect to this 
parameter, they found evidence for a first-order phase 
transition between active (equilibrium fluid) and inactive 
(glass) states. In the second study, Pitard et al. 21 used 
an alternative measure of activity, based on the steep- 
ness and curvature of the energy landscape, integrated 
over time. Using this activity measure to bias the sys- 
tem, they again found evidence for a dynamical phase 
transition, but the properties of the dynamical phases 
were different to those found in Ref. 20, including appar- 
ently non-extensive behaviour of the activity in one of 
the phases. 

In this study, we combine measurements of the differ- 
ent measures of activity used in Ref. 20 and 21. We find 
that these measures couple to different kinds of molecular 



motion. Further, the two measures are anti-correlated in 
the system that we consider. Physically, this happens be- 
cause stable states with very slow structural relaxation 
may have an increased propensity for 'vibrational' mo- 
tion (or /3-relaxation) on short length scales. Based on 
this observation, we are able to resolve some of the appar- 
ent differences between the results of Ref. 20 and 21. We 
also gain insight into the nature of the inactive (glassy) 
states, and how these relate to properties of the underly- 
ing energy landscape, and the normal modes associated 
with motion on this landscape. 

Section II of this paper introduces the model and the 
ensembles that we will use; in Sec. Ill, we compare the 
two measures of the activity used in Ref. 20 and 21, show- 
ing that they are anti-correlated. In Sec. IV, we inves- 
tigate the activity of Pitard et al. 21 in more detail, and 
discuss the relationship of this activity measurement to 
other properties of the fluid and glassy states in the sys- 
tem. We summarise our main conclusions in Sec. V. 



II. BACKGROUND 
A. Model 

We consider the Kob-Andersen mixture of Lennard- 
Jones particles 22,23 , which is a well-studied model glass 
former. There are A particles in the system and a config- 
uration r N has potential energy E(r N ) — J2i<j V( r ij)> 
where is the distance between particles i and j, and 
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There are two species of particle, A (large) and B (small) 
and the parameters and er,;. 
of particles i and j, as ctaa = 

CAB = 0.8(7, £AA = 6 = 1, £BB 

For numerical efficiency, Vij(rij) is truncated at 
2.5<Jij and shifted so that the energy of a pair of particles 



depend on the species 
(7 = 1, <7BB — 0.88cr, 
= 0.5e and 6ab = 1.5e. 

„cut 



separated by r™ is zero. For a system of A particles, 
there are Aa = (4A/5) particles of type A and Ab = 
(A/5) of type B. The density is fixed at p = 1.2(7 -3 as 
in Ref. 20: note that p = 1000/(9.4cr) 3 w 1.204cr -3 was 
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used in Ref. 22 and 23 and in some other studies. This 
small difference has no qualitative effect on the behaviour 
shown here. 

The system evolves by Monte Carlo (MC) dynam- 
ics: as discussed by Berthier and Kob 24 , this dynam- 
ical scheme results in structural relaxation that is in 
quantitative agreement with molecular dynamics, up to 
a rescaling of time. It was also shown in Ref. 20 that MC 
dynamics and constant-temperature molecular dynamics 
gave very similar results in the s-ensemble. The MC dy- 
namical scheme corresponds to a system evolving with 
ovcrdampcd Langcvin dynamics, 



■Vi(t), 



(2) 



where D is the (bare) diffusion constant of a single par- 
ticle, f3 = 1/T is the inverse temperature (we take Boltz- 
mann's constant fee = 1), and rji (t) is white noise with 
zero mean, and covariances 



W{tX{t'))=2D 5 ii 6^5{t-t') 



(3) 



in which fi and v label cartesian components of the vec- 
tor rf(t). The natural units in the system are the length 
a (the diameter of a large particle) ; the energy e (interac- 
tion strength between large particles) ; and the time scale 
At = a 2 /D (of the order of the Brownian time for a free 
particle). When discussing our numerical results in the 
following sections, we take (a, e, At) all equal to unity, 
for compactness. 

The MC dynamical scheme that we use is equivalent 
to the Langevin equation (2) in the limit when all MC 
steps are small (see for example Ref. 25). As in Ref. 24, 
we draw trial MC displacements from a cube of side 
S = 0.15ct, centred on the origin. This choice of step size 
leads to efficient simulations which accurately capture 
the nature of the structural relaxation. The mean square 
displacement for a trial MC move is <5 2 /4: the require- 
ment that the diffusion constant be Do = <J 2 j At means 
that At corresponds to 24 (a/ 5) 2 w 1070 MC sweeps. 

We emphasise that overdamped dynamics as studied 
here were used by Hedges et al. 20 , who also considered 
molecular dynamics with a strong coupling to a thermo- 
stat. However, the results of Pitard et al. 21 were obtained 
using molecular dynamics at constant energy. 



B. Ensembles of trajectories, and measures of activity 

We consider dynamical transitions that occur in en- 
sembles of trajectories. These trajectories have duration 
tobs, and each trajectory is divided into M "slices", each 
of duration At. Following Hedges et al. 20 , the activity of 
a trajectory r N (t) is defined as 



N A M 



K[r N (t)}=AtJ2J2\ r ^- r ^- 



(4) 



where the index i runs over all particles of type A, 
and the tj are the times that separate the slices: tj = 
jAt. We also define the intensive "activity density" 
k = K/(NAt bs), which we sometimes refer to simply 
as the activity. 

From (4), it follows that k measures the mean square 
displacement of a type-A particle during a time interval 
At. This time scale is comparable with the time taken 
for a free particle to diffuse over its own diameter; in the 
supercooled state then At is long enough for a particle 
to explore its local environment (part of the /3-relaxation 
process), but At is shorter than the typical time for the 
fluid structure to relax (the a-process). Our interpreta- 
tion is that k measures motion on length scales compa- 
rable to the particle diameter. 

The dynamical phase transitions that we will consider 
occur when the equilibrium ensemble of trajectories is bi- 
ased to low activity. We define a biased ensemble (or 's- 
ensemble') through its probability distribution over tra- 
jectories: 



P a [r N (t)]<xP [r N (t)]e- K l r "M, 



(5) 



where Po[r N (t)] is the equilibrium probability of trajec- 
tory r N (t). (In defining the probability distributions over 
trajectories, it is sufficient for our purposes to represent a 
trajectory as the set of M + 1 configurations at the times 
tj that separate the slices. However, a finer-grained rep- 
resentation in time is also possible.) 

Within the s-ensemble the average of any trajectory- 
dependent observable A may calculated using 



(c- sK )o 



(6) 



where (-) s denotes an average over trajectories of length 
tobs in the s-ensemble and (-)o means an average of tra- 
jectories of length t b s at equilibrium (which corresponds 
to s = 0). 

An alternative measure of the activity was proposed 
by Pitard et al. 21 , as the time integral (between t = 
and t = i bs) of an 'effective potential': 



y cff = §5>,| 2 + i]>>.F 4 , 



(7) 



where the index i runs over all particles and Fi = —ViE 
is the force on particle i. 
In this study, we define 



At 



M 



i=l j=0 



KMr N (t)] = [^efffo-i) + Kfffe)] (8) 

which is an estimate of the integral of V e s, using a trapez- 
ium rule (we take tj = jAt as above). The notation K a \ t 
indicates that this is an 'alternative' activity. We also 
define fc alt = K &lt /(Nt ohs ), by analogy with k. Since 
V e s is evaluated at only M + 1 points within the trajec- 
tory, K a n is not a very precise estimate of the integral of 
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Kff proposed by Pitard et al. 21 as an activity measure. a ) -350 



However, we expect that K a \ t captures the same physical 
features as this measurement. We also performed simu- 
lations where 2M + 1 points were used to calculate K^t 
(the step size in the trapezium rule was halved). This 
produced no qualitative difference in the values of K a \t 
we obtained. This means that M + 1 points are sufficient 
to make K a \ t a good estimate of the integral of V e g. 

The relation between K a \ t and dynamical activity is 
not obvious a priori. Pitard et al. 21 identified K a \ t as 
an activity by considering the probability that a parti- 
cle returns to (or remains at) its original position over 
a small time St. This probability is obtained from the 
propagator G{r' N ,t';r N ,t) which gives the probability 
that a system in configuration r at time t will evolve 
into configuration r lN at time t' '. For Langevin dynamics 
as considered here, the probability that the initial and fi- 
nal states are the same is given by Autieri et al. 26 : for 
small St, 



G(r N ' ,t + St;r N ,t) = z 



N 



-PV oa 8t+0(St 2 ) 



($ t )3N/2 



= z - 1 e -(3N/2)logSt-l3V o{{ St+0(5t z ) 

(9) 

where z is a normalisation constant (independent of 
time) . We include the full dependence of G on St to em- 
phasise that G decreases with St, regardless of the sign of 
V e fi. (On setting V e g — 0, one recovers the standard re- 
sult for N non-interacting Brownian particles.) Equ. (9) 
shows that when V e g is large then particles in the sys- 
tem are likely to move quickly away from their original 
positions; when V c g is small then particle are more likely 
to remain localised. This is the motivation for proposing 
K a \ t as a measure of dynamical activity. Note however 
that this measurement is defined in terms of motion on 
the very small time scale St. 

Following Pitard et al? 1 , we therefore define an 's a i t - 
ensemble' through a bias on K a \ t : 

P Salt [r N (t)} oc P [r N (t)}e- s ^ K ^ [rN{t)] . (10) 

This definition is analogous to (5): continuing the anal- 
ogy for averages of an observable A, we have 



(11) 



by analogy with (6). Equations (5) and (10) define the 
ensembles of trajectories that we will consider in the fol- 
lowing. 



III. MEASUREMENTS OF ACTIVITIES IN BIASED 
ENSEMBLES 

We use transition path sampling (TPS) 28 to sample 
biased ensembles of trajectories, as discussed in Ap- 
pendix A. We show numerical results obtained by TPS in 
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FIG. 1. (a) Scatter plot of the two activity measurements k 
and fc a it, in three different s-ensembles. The ensembles are 
characteristic of the active phase (s — 0.000), the coexistence 
region (s = 0.015) and the inactive phase (s — 0.025). The 
two activity measurements k and fc a i t are anti-correlated. The 
trajectory length is t b s = 400At. (b, c) Marginal distribu- 
tions of k and fc a i t from the s-ensemble with s — 0.015. This 
bimodal behaviour is characteristic of the dynamical phase 
transition found in Ref. 20. (d) Scatter plot of k and fc a it 
for three values of s a it and f b s = 200At. The data for 
Salt = —3.0 x 10~ 5 is similar to the inactive data for s — 0.025. 
The dashed and dotted lines in (a) and (d) are the same in 
both panels and are obtained by linear regression analyses on 
data from (a) for the dots and (d) for the dashes. 



Figs. 1 and 2, which summarise the behaviour of K and 
K a \ t , as s and s a it are varied. We concentrate on the be- 
haviour of a system of N = 150 particles at temperature 
T = 0.6, as in Ref. 20. [Recall we have fixed units such 
that (e,a,At) are all equal to unity] In Figs. l(a,d), we 
show scatter plots of K and K a \t, combining data sam- 
pled from equilibrium and for several values of s and Salt- 
We find that k is always positive and fc a it is always neg- 
ative. (It is worth noting that throughout this article we 
write "fc a it is larger than k' alt " if fc a it > fc a it> regardless 
of the sign of fc a it-) Perhaps surprisingly, we also find 
that while k and k a \ t were both proposed as measures 
of dynamical activity, they are anti-correlated with one 
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FIG. 2. Averaged activities in biased ensembles. Note 
that panel (b) shows the negatives of the field and the ac- 
tivity, — s a it and (— feait). All data are for N = 150 and 
T = 0.6, except for the red dashed lines, where N = 300 
and we show the linear response behaviour about equilib- 
rium: (K) s = (K)o + s(6K 2 )o + 0(s 2 ), and similarly for s a it- 
These linear response results do not capture the non-trivial 
crossovers, but they do show that the mean and variance of K 
and Ka\t are approximately extensive in N, for s = (there is 
a weak finite-size correction to (fc)o: particle motion in smaller 
systems is known to be slightly slower for this system, com- 
pared to the bulk 27 ). 

another. This observation will be crucial in the following 
discussion. 

Panels (b) and (c) of Fig. 1 also show that for an ap- 
propriate value of s (here s — s* = 0.015), the marginal 
distributions of both fc and fc a it are bimodal. These dis- 
tributions are indicative of a dynamical phase transition, 
although the existence of such a transition can be con- 
firmed only if these distributions remain bimodal as the 
system size TV and observation time t Q b s are taken to in- 
finity. 

In Fig. 2, we show average values of k and fc a it in 
ensembles of trajectories, as s and s a it are varied. We 
note however that in Fig. 2(b), we are plotting (— fc a it) 
against — s a it . On increasing s in panel (a), we observe 
a crossover from a large-fc state at s = to a small- fc 
state at positive s. As in Ref. 20, this crossover becomes 
sharper as i bs is increased, consistent with a dynamical 
first-order phase transition. As we increase —Salt (or de- 



crease Salt) in Fig. 2(b), we observe a similar crossover to 
a state with smaller — fc a lt (and hence larger fc a it)- Again, 
the crossover sharpens on increasing i Q b s - 

Finally, returning to Fig. l(a,d), we observe that the 
states for s = 0.025 and s a it = —3 x 10~ 5 have simi- 
lar joint distributions of (fc,fc a it)- Hence, taking Figs. 1 
and 2 together, we infer that the two crossovers shown in 
Fig. 2 represent transitions between the same two states: 
the equilibrium state [colored red in Fig. l(a,d)] and the 
state that was identified by Hedges et al. as the glassy (in- 
active) state [colored blue in Fig. l(a,d)]. It was shown by 
Hedges et al. 20 that the inactive state was accompanied 
by a self-intermediate scattering function that does not 
decay throughout the observation time i bs, indicating 
that particles remain localised near their initial positions 
throughout the trajectory. Our data confirm this result: 
this is the sense in which this small-fc state is 'inactive'. 

The crossover shown in Fig. 2(b) for small negative s a it 
was not reported by Pitard et al. . However, we note 
that the ranges of s a it and K a \ t shown in Fig. 2(b) are 
much smaller than those used in Ref. 21. It is possible 
that a more detailed analysis of the relevant range of s a it 
using the methodology of Ref. 21 might reveal a similar 
crossover/transition. What is clear from Figs. 1 and 2 is 
that the transition (for s a it > 0) reported by Pitard et 
al 21 is a different phenomenon to that reported by Hedges 
et al 20 . 

The transition reported by Pitard et al 21 for s a it > 
is accompanied by anomalous behaviour of the deriva- 
tive dfc a it / ds a it and non-extensivity of fc a it itself, for small 
positive Salt (and perhaps even for s a it = 0). For the nar- 
row range of s a it that we considered, we did not observe 
these effects. Fig. 2(b) shows that dfc a it/ds a i t depends 
very weakly on s a it for s a it < 0, and that on increasing 
the system size to 300 particles, there is no significant in 
change either the equilibrium average of fc a it nor in its 
derivative with respect to s a it- The differences between 
our results and those of Ref. 21 in this regime remain 
a subject for future study: here we concentrate on the 
crossover that we do find for s a it < 0, and its relation- 
ship to the active/inactive phase coexistence phenomena 
found in Ref. 20. 



IV. INTERPRETATION OF ACTIVITY 
MEASUREMENTS 

The interpretation of the activity fc is transparent in 
that it measures particle motion on a timescale At. As 
discussed in Ref. 20, the low-fc phase found on increasing 
s is characterised by an absence of structural relaxation 
(at least for small systems of 150 particles, on time scales 
up to 40 times the equilibrium relaxation time). The 
relation between fc a i t and particle motion is somewhat 
indirect, operating via the expression (9) which gives the 
probability that a particle deviates significantly from its 
initial position, on short time scales. 

In the following, we focus on the activity fc a it, aiming in 
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FIG. 3. Numerical test of equation (12). The two quanti- 
ties should be equal at equilibrium (s — 0), but there is a 
small difference between them due to our use of a truncated 
potential (see Appendix 12). The small difference is almost 
constant for the range of s considered. 

particular to understand why this activity measurement 
is larger in the 'inactive state' of Ref. 20, compared with 
equilibrium. 



A. Two contributions to V c s, and a 
quasi-equilibrium/two-temperature scenario 

From (7), we see that V e g (and hence also k a \ t ) has two 
contributions, one from the interparticle forces and the 
other from the divergence of the force. At equilibrium, 
these contributions are related: 

(m\\ =Z- 1 j dr N \l3S7 t E(r N )\ 2 e^ E ^ 

=Z~ l J dr N (3V 2 l E(r N )e^ E ^ {l2 ^ 
= - ((3V, ■ F t } 

where Z = J dr N c~ /3E ( rN ' > is the equilibrium partition 
function. The first and third equalities in (12) fol- 
low trivially from the definition of the equilibrium av- 
erage, while the second relies on an integral by parts. 
This result is well known and has been exploited to de- 
termine the temperature of a system directly from its 
configurations 29 ' 30 . At equilibrium, we conclude that 

(K ff )o--!E«(l^l 2 )o. 

Data for the two terms in V e s is shown in Fig. 3. De- 
spite (12), we note a small difference between the two 
terms, even at equilibrium. This effect arises because of 
the truncated and shifted Lennard- Jones potential that 
we use in simulation, which has a discontinuity in its first 
derivative at the cutoff radius r™* = 2.5<7ij. 

We discuss this effect in Appendix B (see also 30 ) where 
we define a regularised average (V • Fi) slm , and discuss 
how (12) is modified to account for this regularisation. 
Consistent with Fig. 3, we find that the effect of this 



regularisation is small throughout, so we use (V • i*i) slm 
interchangeably with (V • Fj) in what follows. 

Having accounted for the small systematic deviation 
between the two quantities plotted in (3), the most im- 
portant feature of that figure is that the two contribu- 
tions to V e fi remain almost equal, as s increases. That 
is, for the range of s considered, our numerical results 
indicate that 

fait), « -^(|F 4 | 2 ) S « l^V* (13) 

i i 

Since (12) applies only at equilibrium, this is a non-trivial 
result. Our interpretation is that the 'slow' (structural) 
degrees of freedom respond strongly to the bias s, while 
the 'fast' (or 'vibrational') degrees of freedom respond 
much more weakly. In other words biasing moves the 
system to a region of the energy landscape not typical of 
equilibrium, but the system explores that region as if it 
were at equilibrium. If this is indeed the case, the equilib- 
rium assumption required to prove (12) can be replaced 
by a weaker, 'quasi-equilibrium' assumption for the fast 
modes, leading to a similar result. 

We formalise this hypothesis within a mean-field 
description 5 ' 31 , assuming that the system has many 
metastable states. A short relaxation time is associ- 
ated with intrastate ( "vibrational" ) motion and a longer 
relaxation time is associated with structural rearrange- 
ment (between states) 10 . We emphasise that metastable 
states are defined dynamically, by reference to their life- 
time 10 ' 32 : each state contains many energy minima ('in- 
herent structures' 7 ). 

Following the discussion of Ref. 10, for a weak bias s 
then the steady state distribution over configurations C 
is 

p ss (C)^w a(c)C -f 3E ^/Z a(c) (14) 

where a(C) is the state containing configuration C, 
while w a (c) is the probability of that state, and Z a = 
Scea e~^ E ^ is the equilibrium weight of state a. (Here, 
the short-hand notation C indicates a configuration r N .) 
The s-dependence of (14) comes only from the weights 
w a . If w a = Z a for all states a then we recover the 
equilibrium Boltzmann distribution (at s = 0). For fi- 
nite s then one expects the w a associated with long-lived 
metastable states to be enhanced. A similar idea was 
discussed in Ref. 33, where inherent structures were used 
in place of metastable states. 

As usual with mean-field scenarios, Equ. (14) is ap- 
proximate for (at least) two reasons: firstly, it assumes 
that each configuration can be assigned to a single 
metastable state (which neglects configurations on the 
boundaries between states); secondly it assumes that 
intra-state fluctuations are unaffected by the field s. The 
first approximation can be ignored in mean-field mod- 
els because configurations on boundaries between states 
have negligible weight in p ss (C). The second approxima- 
tion is valid for small s, if (and only if) fast and slow 
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dynamics take place on well-separated time scales. This 
situation is realised in mean-field models and may be ex- 
pressed in terms of a condition on the eigenvalues of the 
time evolution operator of the system 10 . 

In finite-dimensional systems (where mean-field theory 
is not exact), both of these approximations lead to devi- 
ations from (14), but one expects that equation to give a 
reasonable description of the system if the lifetimes of the 
metastable (inactive) states are much longer than time 
scales for motion within these states. Ref. 33 shows that 
this condition is quite well-satisfied. Hence, one may re- 
peat the analysis of Equ. (12), but using (14) in place 
of the Boltzmann distribution. One arrives at the same 
conclusion, that the two terms plotted in Fig. 3 should 
be equal. The largest error in that analysis comes from 
configurations that lie on boundaries between metastable 
states 34 , but our numerical results indicate that these 
configurations do not contribute too much to these aver- 
ages, and that the quasi-equilibrium hypothesis of (14) 
seems to hold quite accurately. This is the sense in which 
the slow fluctuations (between states) respond strongly 
to the field s (via the w a ), while the fast (intra-state) 
fluctuations respond much more weakly. 



B. Vibrational modes of the fluid in biased ensembles 

The relationship between k a n and the properties of 
equilibrium and inactive states can also be analysed 
through the distribution of eigenvalues of the dynamical 
matrix (or Hessian) H . This distribution, together with 
the vibrational normal modes of supercooled liquids have 
been connected with their dynamical properties in a va- 
riety of studies 17-19,35 . Here we exploit the connection 
between the matrix H and the contribution of Vj • Fi 
to fc a it. The Hessian is a 3N x 3N matrix with elements 
Hi^i.jv = ffygfy i where the indices i and j run over all 
particles and fi and v run over the cartesian components 
of the position vectors r». 

The matrix H has 3N eigenvalues, which we denote by 
oj 2 ,u>2, Here, each uj a can be interpreted as a natu- 
ral frequency for vibrational motion on the energy land- 
scape, along a particular eigenvector. However, we note 
that since typical configurations of the system are not 
located at minima of the energy landscape, some eigen- 
values of H will be negative, io 2 a < 0. In this case the 
interpretation of w a is less clear, but the relevant direc- 
tions on the energy landscape are unstable, indicating 
that the system is close to a saddle point of the land- 
scape, and not a stable minimum. The V • F term in V e g 
is related to the eigenvalues as 

3N 

£v i .fl = -'ft(S) = -ZX (15) 

i a—1 

Defining the distribution of eigenvalues, D(uj 2 ), the 
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FIG. 4. (a) Distribution of eigenvalues of the Hessian for both 
phases, (a, inset) The difference AD(cj 2 ) = [D(u; 2 ) s= o.o4 — 
D(ui 2 ) s =o.oo] between the phases. The distribution for the 
active phase is slightly broader, and it associated mean value 
of cu 2 is larger, (b) Distribution of oj where cu 2 > for both 
phases. 

trace can be expressed as 

/oo 
d{Lu 2 )uj 2 D(u 2 ) (16) 
-oo 

Combining (13) and (15) and (16), we see that k a it ~ 
~~T^ f-oo^ 1 ^ 2 ) w2 -D(w 2 ), allowing us to relate the differ- 
ence in k a it between active and inactive (small-fc) states 
to the distribution D(ui 2 ) of these states. Results are 
shown in Fig. 4. Comparing equilibrium (s = 0) and in- 
active (s > 0) data, the differences in D(ui 2 ) are subtle, 
but the dominant effect is that the main peak in D{lo 2 ) 
is slightly sharper in the inactive state. That is, the in- 
active state has fewer modes with small or negative u> 2 , 
but also fewer modes with large positive cu 2 . Hence it 
has more modes with intermediate lo 2 . When evaluat- 
ing the change in Tr(H) between states, the dominant 
effect comes from large eigenvalues, which correspond to 
"stiff" (strongly-curving) directions on the potential en- 
ergy landscape. Fig. 4 shows that there are fewer stiff 
directions in the inactive state, and this results in k a \t 
being larger (less negative) for that state. The difference 
is more pronounced when plotting _Di(w), the distribu- 
tion of to among modes where lo 2 > 0. 
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FIG. 5. (a) Distribution of eigenvalues of the Hessian for 
inherent structures of both phases, (b) Distribution of uj for 
inherent structures of both phases, (b, inset) Dividing Di(u>) 
by uj 2 emphasises the lack of low frequency modes associated 
with the inactive phase. 
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FIG. 6. A schematic representation of the differences in 
the energy landscape between the active and inactive phase. 
In the inactive phase, the barriers between basins (inherent 
structures) are smaller making rearrangements on large length 
scales less likely. These correspond to small values of uj 2 . The 
strongly curving directions around basins are less steep in the 
inactive phase, allowing more motion on short length scales. 
These correspond to large values of ui 2 . 



In Fig. 5, we show the distributions of uj 2 and of uj that 
we obtained by using conjugate gradient minimisation on 
configurations from the s-ensemble, and then construct- 
ing the matrix H at the resulting energy minimum [in- 
herent structure (IS)]. In this case, all eigenvalues of H 
are positive. The differences in D(uj 2 ) between active and 
inactive states are more pronounced at the IS level, but 
the main conclusion is the same: the peak in D{uj 2 ) is 
narrower in the inactive state, and this pushes the mean 
value of uj 2 to a smaller value. However, these data also 
emphasise that the inactive state has fewer "soft" modes 
(with small uj), compared to equilibrium. This effect was 
noted in Ref. 33: it indicates that part of the stability of 
the inactive state can be accounted for by the paucity of 
soft-directions on the energy landscape. 

The resulting physical picture is summarised in Fig. 6. 
The potential energy surface (or 'landscape') is divided 
into basins, each associated with a single inherent struc- 
ture (local minimum). Moving away from the inherent 
structure, most of the directions are quite 'stiff', with 
large uj, but a few are 'soft', with small uj. Comparing 
the equilibrium state with the inactive (small-fc) state, 
Figs. 4 and 5 show that the stiff directions in the inac- 
tive state are (on average) less stiff than at equilibrium; 
on the other hand, the soft directions in the inactive state 
are also less soft than at equilibrium. The activity pa- 
rameter fe a it of Pitard et al. 21 is most sensitive to the stiff 
directions: the stiffer these are, the less particles are free 
to move (on short scales), and the smaller is /c a it- On the 
other hand, the activity parameter k of Hedges et al. 20 
is most sensitive to structural relaxation, which couples 
more strongly to the soft modes: these are less soft in 
the inactive state, suppressing large-scale particle mo- 
tion, and reducing k. 

This difference in sensitivity to fast and slow motion 
explains the anticorrelation between k and fc a i t in Fig. 1, 
and it also explains why the active/inactive transition 
of Ref. 20 appears only in s a it-ensembles with s a it < 0. 
We argue that it should be borne in mind in any future 
studies that use V e ff to measure activity. 



C. Liquid structure in biased ensembles 

We now turn to the structure of the active and inactive 
states that we have found, and the connection of this 
structure to fc a it- It is notable from Fig. 2 that typical 
values of fc a it are around — 380(e/er 2 ), while the difference 
in fc a it between active and inactive states is much smaller, 
around 30 (e/er 2 ). (We give the units of fc a it explicitly 
in this discussion: recall that numerical data are shown 
after fixing (e, a) to unity.) 

To interpret these results, it is useful to write 



( E V ' F *)s = E / 4 7 rr 2 dr5 ij -(r)V 2 ^(r) 



(17) 



where gij(r) = (J2j^i^( r ~ r ij))s is proportional to a 
radial distribution function (in the s-ensemble). Since 
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Vij(r) and (jij(r) depend on the particle indices I and j 
only through their types, it is convenient to use a short- 
hand notation for the non-trivial part of the integrand in 
(17) 



a) 



G (r) = V 2 Vij {r)(jij (r) 



i,j of type A 



(18) 



where the right hand side is evaluated with i and j both 
being particles of type A. Similarly, we define G AB (r) 
and G BB (r) for particles of other types. (Note that these 
functions depend implicitly on the biasing parameter s, 
through^-.) 

By comparing 4nr 2 G AA (r) to <7aa(?") (the radial distri- 
bution function for particles of species A) , we can see how 
the liquid structure on different length scales contributes 
to (V.Fi) s . We focus only on the function for the large 
particles as these are the most numerous species. 

Fig. 7 (a) shows saaM for the active phase (at 
s = 0.00) and the inactive phase (at s = 0.04). There are 
some subtle changes: the first and second peaks and the 
first trough are enhanced in the inactive phase. Panel 
(b) shows 47rr 2 G AA (r) for the same values of s. Only 
a small region contributes to (V.-Fj) - the width is less 
than that of the first peak in Qaa{t). This further em- 
phasises that fc a n is dominated by behaviour on short 
length scales. Again, the differences between the phases 
are subtle. This is in line with the observation that the 
size of the change in fc a i t between phases is much smaller 
than the size of k a \t itself. 

To emphasise the change, we consider the difference 
AG AA (r) = [G AA (r)] s=0 . 04 - [G AA (r)] s=0 . 00 . This is 
shown in the inset to figure 7 (b). It is clear that the 
change in fc a i t is largely due to changes in the liquid 
structure at very small length scales; the dashed line in 
the plot indicates where G AA (r) is largest in magnitude, 
which corresponds to the maximum of the first peak in 
S'aa(^)- These changes are subtle enough that they are 
not apparent when comparing radial distribution func- 
tions, but since V 2 I4,(r) is very large for small r they 
are ultimately what is important when considering fc a i t . 

In addition to the results in Fig. 7, we have obtained 
similar data for G AB (r) and G BB {r): the main picture 
is the same but the smaller numbers of B particles in the 
system mean that these functions contribute less strongly 
to V e g, and also that the numerical uncertainties in our 
results are larger. As shown by Speck and coworkers 36,37 , 
the radial distribution function g BB (r) shows the largest 
relative changes between active and inactive states. How- 
ever, the small number of B-particles means that this 
gives a relatively small contribution to the changes in 
fc a it shown in Fig. 2. 



D. The dynamical action 

Finally, we discuss one other context in which the ac- 
tivity -ftT a it appears. For overdamped dynamics as in (2), 
at equilibrium, the probability of a trajectory r N (t) can 
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FIG. 7. (a) Comparison of the partial pair correlation func- 
tion for large particles between the active and inactive phase. 
Although there are some differences (the height of the first 
peak and the depth of the first trough) they are small, (b) 
The function 4vrr 2 G AA (r) which can be integrated to give 
(V • Fi) s . The interesting part of this function occurs around 
the position of the first peak in the pair correlation function. 
The inset panel shows the difference in this function between 
the phases, AG AA (r) = [G AA (r)] s=0 . 04 - [G AA (r)] 8=0M) . This 
serves to illustrate that the changes in A" a it come from struc- 
tural changes on short length scales. 



be written as 26 



Po[r"(t)] = -P il ee[r"(t)]-et 



\[E(0)-E(t ob .)] 



-fSD K^ t [v N {t)} 



(19) 



where Pf rcc [r N (t)] is the probability of the trajectory in 
the absence of any forces, and Z is a normalisation con- 
stant. 

Hence if we consider the equilibrium distribution of fc a i t 
for this model, we have 

iWfc.lt) = le^ 5 ^)-^-], (20) 

where e JV * obs<s ( fcalt ) is the marginal distribution of k & \ t as- 
sociated with the distribution P bcc [r N (t)]e^ [Em ' E ^° b ^ ] . 
(We emphasise that the function tS(fc a it) depends on the 
parameter f3 via the definition of fc a it , and it also depends 
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on D .) Further, the distribution of k a \ t within the s a i t - 
ensemble is 

P s (fcait) oc e m ^ s ^-^ Do+s ^ k ^\ (21) 

There is a relevant analogy here: compare the distribu- 
tion of the energy density e = E/N in a thermal system 
at equilibrium, 

P fj (e)o,e N ^-M, (22) 

where S(e) is the entropy per particle. This analogy be- 
tween ensembles of trajectories like (21) and ensembles 
of configurations like (21) was a key starting point for 
studies of the dynamical transitions and biased ensem- 
bles that we consider here 9 ' 12-14 . 

Extending this analogy, the interpretation of k a i t and 
Salt is as follows. Within the distribution Po[r N (t)], 
there are many trajectories with large values of fc a it, 
each of which is individually rare because of the factor 
of e~P D o K ^. There are fewer trajectories with smaller 
^aitj but these are individually more probable because 
they are less strongly suppressed by the factor e~ l3D " ka - lt . 
The most likely value of fc a it occurs when the 'entropic' 
term iS(fc a it) balances the 'energetic' term f3Dok a \ t . [Here 
we are using the labels 'entropic'/'energetic' to empha- 
sise the analogy with (22): these terms have no simple 
relation to thermodynamic energy or entropy] 

If we introduce a negative value of s a i t , the sys- 
tem is biased towards the more numerous ('entropically 
favourable') trajectories in the system, which have larger 
(or less negative) values of fc a i t . As shown in Fig. 2, even 
a small negative s a it is sufficient to drive the system into 
an 'inactive' state in which structural relaxation is ar- 
rested. The unexpected anticorrelation between k and 
fc a it that we found in this study arises because the inac- 
tive state has the higher 'entropy' S in trajectory space. 
The reason for this is that the inactive state consists of 
configurations in in which most directions on the energy 
landscape are not too 'stiff': despite the slow structural 
relaxation, the particles have greater freedom to move on 
small length scales, compared with equilibrium. And the 
more free the particles are to move, the more trajectories 
are available, and the larger is S. As before, the con- 
clusion is that propensity for motion on small scales is 
anti-correlated with propensity on scales of the order of 
the particle diameter. 

V. CONCLUSIONS AND OUTLOOK 

This study has two central conclusions. Firstly, the 
transition found by Hedges et al. 20 for s > corresponds 
to a transition for s a \ t < within the ensembles defined 
by Pitard et al. 21 . Secondly, the activity parameter k a it 
defined in Rcf. 21 couples to dynamical motion on small 
scales, which is anticorrelated with the structural relax- 
ation of the fluid. This anticorrelation arises from prop- 
erties of the energy landscape of the inactive state. In 



addition to these main points, we have also discussed the 
structure of the inactive states and the connection of k a \ t 
the liquid structure; and also the extent to which the in- 
active states have the quasi-cquilibrium property given 
in (14). 

We hope that this work clarifies the role of the activ- 
ity measurement introduced by Pitard et al. 21 , which we 
have denoted by K a \ t . Equ. (21) shows that K a \ t is inti- 
mately connected with dynamical motion in overdamped 
Langevin systems, and it is also strongly connected to 
the energy landscape of the fluid. These facts present a 
strong argument in favour of K a \ t as an activity measure 
that arises naturally from the dynamics of the system, 
without any prejudice as to the nature of its dynamical 
relaxation. However, the results of Fig. 1 show that K a \ t 
must be interpreted carefully, since the extent of short- 
scale motion may not be correlated with the effectiveness 
of structural relaxation. Also, this study did not find ev- 
idence for singular behaviour in (fc a it) Salt for the range of 
positive s a it that we considered: the physical interpreta- 
tion of the behaviour found in Ref. 21 for larger positive 
s a it remains unexplained (although it seems unrelated to 
the active/inactive crossover discussed in Ref. 20). 
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Appendix A: Sampling biased ensembles 

We sample trajectories from the s-ensemblc and s a i t - 
cnscmblc by using transition path sampling (TPS). This 
method samples trajectories in a similar way to the sam- 
pling of configurations by standard Metropolis Monte 
Carlo methods. Its operation is reviewed in Ref. 28 and 
the 'shifting moves' used in this study are discussed in 
Ref. 38. We give a brief overview here: Starting with 
an initial trajectory r$ (i), a new trajectory rf(t) is 
generated by a 'shifting move'. In 'forward shifting', 
one chooses a random number p between 1 and M, 
and slices 1,2, ...,p of r^(i) are discarded. The re- 
maining slices (p + 1 , . . . , M) of (t) form the initial 
slices (1, . . . , M — p) of the new trajectory rf (t). Slices 
M — p + 1 , . . . , M are then generated by unbiased dy- 
namical evolution from slice M — p. Finally, this new 
trajectory rf (t) is accepted with probability 

Pace - min{l,e-^[^ t )]+^^ t )]} . (Al) 

Otherwise one rejects the new trajectory and retains the 
original one, (t). This procedure is used in conjunction 
with "backwards shifting" moves where slices 1,2, ... ,p 
of (t) are used to form slices M — p + 1,...,M of 
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r^(i), and then slices 1, . . . , M — p of r^(t) are gener- 
ated by unbiased time evolution, backwards in time from 
slice M — p + 1 (use of this scheme requires the time- 
reversal symmetry property of the equilibrium state of 
the model) . This combination of moves ensures detailed 
balance within the ensemble of trajectories (5), so after 
sufficiently many moves, the procedure converges in a 
stationary regime which generates representative samples 
of the ensemble. Further, since the system is stochastic 
and the ensemble of trajectories being sampled is (ap- 
proximately) time-translationally invariant, these shift- 
ing moves are effective in sampling the ensemble, and it is 
not necessary to supplement them with 'shooting' moves. 
(A combination of shooting and shifting is the conven- 
tional choice in rare event sampling problems dominated 
by barrier crossing, but we do not use this procedure 
here). 

The results shown here were obtained from TPS simu- 
lations as follows. We used a weighted histogram analysis 
(WHAM) 39 to combine data obtained using different val- 
ues of s and s a \ t . For trajectories of length t Q b s = 200 At 
we used data from s = —0.025 to s = 0.03 in the s- 
ensemble and from s a \ t = —3.0 x 10~ 5 to s a i t = 5.0 x 10 -5 
in the s a i t -ensemble. For trajectories of length t b s = 
400A£ we used data from s = 0.00 to s = 0.020 for the s- 
ensemble and from s a \t = —1.75 x 10~ 5 to s a i t = 0.00 for 
the s a it-ensemble. These choices ensure that we concen- 
trate our numerical effort in the crossover regime between 
active and inactive states: as we bias further into the in- 
active regime, the slow structural dynamics of the inac- 
tive state limit the effectiveness of sampling. We there- 
fore access the inactive regime by histogram reweighting 
from the crossover regime, using the results from WHAM. 

Large values of s (and — s a it) bias the system towards 
inactive states, and this can lead to crystallisation within 
trajectories. This happens rarely and we exclude trajec- 
tories with a high degree of crystalline order from our 
analysis. We measure crystalline order using the common 
neighbour analysis scheme described in the supplement 
to Ref. 20. We note that the values given for the max- 
imum separation of bonded pairs of particles in Ref. 20 
are incorrect, and we use the correct values: Aaa = 1-45, 
Aab = 1-25 and Abb = 1-07. 

We note that Pitard et al. 21 used a different method 40 
to sample biased ensembles of trajectories. In contrast to 
transition path sampling, which operates on trajectories 
of fixed duration i bs, that method provides direct esti- 
mates of observables in the limit where t b s oo. On 
the other hand, the algorithm requires that many copies 
(or clones) of the system evolve in parallel, and there are 
systematic errors associated with the method 40 , which 
vanish only when the number of clones is taken to infin- 
ity. In this sense, the TPS method results in controlled 
sampling of ensembles with finite t b s , requiring an ex- 
trapolation to reach the large-i b s limit; on the other 
hand, the method of Ref. 40 gives direct access to a limit 
of large i b s , but at the expense of an extrapolation in 
the number of clones. 



Appendix B: Regularisation of V ■ F> 

The results in Fig. 3 indicate that Eq. (12) is not sat- 
isfied exactly at equilibrium, for the model system used 
here. As discussed in Ref. 30, this behaviour is generic for 
systems where interaction potentials are truncated. To 
analyse this behaviour quantitatively, we imagine modi- 
fying the potential Vij (r^ ) in a region of width e around 
r™* so that its second derivative exists everywhere, and 
then taking the limit of small e. In this case, 

V, • F, = [QH + QiAm rffj] (Bl) 

where 

/ -V%(r«), r y <rff (m) 
13 [0 otherwise, 

and qij = — — is the discontinuity in the force at 
the potential cutoff. If one uses (Bl) as the definition of 
V • Fi, then (12) will hold exactly at equilibrium. 

However, the 6- function in (Bl) makes it problematic 
in simulation. We therefore define instead 

<E v -^) sim = (E^> ( B3 ) 

i i^tj 

and note that 

( ]T V • F t ) = < ]T V • F,) sim + A APA A AA 

i i 

+ AapbAab + AapaAbb (B4) 

where (J^i V-.Fi) on the left hand side uses the definition 
from (Bl), while A AA = 47r(r A u ^) 2 g AA 5 AA (r A u ^),with 
similar expressions for Aab, Abb- Here p\ = Na/V 
is the number density of A-particles, g AA {r) is the ra- 
dial distribution function between A particles, and g AA 
is the value of qij if particles i, j are both of type A. (We 
used the fact that if particles i and j are of type A then 
(5(r — rij)) = 4:Trr 2 pAg AA (r)). We have evaluated the A- 
terms in (B4) at equilibrium, and verified that the data 
in Fig. 3 are then consistent with (12). However, since 
these A- terms are small, we use V-Fi) slm throughout 
this work as our numerical estimator for V • Fi). 
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